### "Elitist Remedies: Complaint Resources and Representation in International Human Rights Bodies" ###
### Replication files for study 1 (national level) ###
### Author: Christoph Valentin Steinert ###

# Packages
library(stargazer)
library(lmtest)
library(pscl)
library(sandwich)
library(car)
library(sjPlot)
library(ggthemes)  
library(texreg)
library(tidyverse)
library(Hmisc)
library(foreign)
library(MASS)
library(lme4)
library(giscoR)
library(sf)
library(countrycode)
library(modelsummary)

# Empty environment
rm(list = ls())

## Load data
load("analysisdata_study1.RData")

# Joint distribution of Count of Communications and Political Terror Scale
### REPRODUCES FIGURE 3 IN THE ARTICLE ###
#pdf(file = "Figure3.pdf")
#tiff("Figure3.tiff", units="in", width=5, height=5, res=1200)
ggplot(pts, aes(PTS, comm_count)) + 
  geom_jitter() + xlab("Political Terror Scale") +
  ylab("UNSP Communications count") +
  geom_smooth(method = "loess", span = 0.5) +
  theme_bw()
#dev.off()
#dev.off()
ggplot(pts, aes(PTS, comm_count)) + 
  geom_point() + xlab("Political Terror Scale") +
  ylab("Communication count") +
  geom_smooth(method = "loess", span = 0.5) +
  theme_bw()
ggplot(pts, aes(PTS, comm_count)) + 
  xlab("Political Terror Scale") +
  ylab("Communication count") +
  geom_smooth(method = "loess", span = 0.5) +
  theme_bw()
### Label outliers
ggplot(pts, aes(PTS, comm_count)) + 
  geom_jitter() + xlab("Political Terror Scale") +
  ylab("Communication count") +
  geom_smooth(method = "loess", span = 0.5) +
  theme_bw() + geom_text(aes(label=ifelse(comm_count>50,as.character(country_year),'')),hjust=1.08,vjust= 0.2)

# Conditional means in PTS categories
tapply(pts$comm_count, INDEX=list(pts$PTS), FUN=mean)
tapply(pts$comm_binary, INDEX=list(pts$PTS), FUN=mean)

# DV histogram and density
hist(pts$comm_count, freq = T, breaks = 50, main = "", xlab = "# of Communications on physical integrity rights abuses")
plot(density(pts$comm_count), main = "")


########################### WORLD MAPS ###########################################

#### Map PTS-scores averaged over 10 years ####
pts$PTS_A <- as.numeric(pts$PTS_A)
pts_averaged <- pts %>% 
  group_by(country) %>% 
  summarise(pts_avg = mean(PTS_A, na.rm = T)) %>%
  filter(pts_avg != "NaN") %>% 
  mutate(pts_rank = rank(desc(pts_avg), ties.method = "min")) %>% 
  arrange(desc(pts_avg)) %>% 
  ungroup()

# Mapping PTS
countries <- gisco_get_countries(year = "2016") # does only work for this year
countries <- subset(countries, CNTR_NAME != "Antarctica")
countries_small <- dplyr::select(countries, ISO3_CODE, geometry)
pts_averaged$ISO3_CODE <- countrycode(pts_averaged$country, origin = "country.name", destination = "iso3c")
pts_averaged <- merge(pts_averaged, countries_small, by = "ISO3_CODE", all.x = T, all.y = T) 
### REPRODUCES FIGURE 2 IN THE ARTICLE ###
ggplot(data = pts_averaged$geometry) + 
  geom_sf(aes(fill = pts_averaged$pts_avg), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "Average PTS-score \n (2010-2020)") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
#ggsave("Figure2.pdf", device = cairo_pdf(), width = 297, height = 210, units = "mm")

# Count number of executions by country and rank
country_based_executions <- pts %>% 
  group_by(country) %>% 
  summarise(execution_count = sum(execution_count)) %>% 
  arrange(desc(execution_count)) %>%
  mutate(execution_rank = rank(desc(execution_count), ties.method = "min")) %>%
  ungroup()
# Count number of detentions by country and rank
country_based_detentions <- pts %>% 
  group_by(country) %>% 
  summarise(detention_count = sum(detention_count)) %>% 
  arrange(desc(detention_count)) %>%
  mutate(detention_rank = rank(desc(detention_count), ties.method = "min")) %>%
  ungroup()
# Count number of disappearances by country and rank
country_based_disappearances <- pts %>% 
  group_by(country) %>% 
  summarise(disappearance_count = sum(disappearance_count)) %>% 
  arrange(desc(disappearance_count)) %>%
  mutate(disappearance_rank = rank(desc(disappearance_count), ties.method = "min")) %>% 
  ungroup()
# Count number of torture by country and rank
country_based_torture <- pts %>% 
  group_by(country) %>% 
  summarise(torture_count = sum(torture_count)) %>% 
  arrange(desc(torture_count)) %>%
  mutate(torture_rank = rank(desc(torture_count), ties.method = "min")) %>% 
  ungroup()

# Merge to one dataset
datacomp <- merge(pts_averaged, country_based_detentions, by = "country")
datacomp <- merge(datacomp, country_based_disappearances, by = "country")
datacomp <- merge(datacomp, country_based_torture, by = "country")
datacomp <- merge(datacomp, country_based_executions, by = "country")

# Count all communications
datacomp$all_commcount <- (datacomp$detention_count + datacomp$torture_count +
                             datacomp$disappearance_count + datacomp$execution_count)
datacomp$allcomm_rank <- rank(desc(datacomp$all_commcount), ties.method = "min")

# Add all countries included in PTS averaged
notincluded <- anti_join(pts_averaged, datacomp, by = "ISO3_CODE")

# Add all missing variables for row binding
notincluded$detention_count <- NA
notincluded$detention_rank <- NA
notincluded$disappearance_count <- NA
notincluded$disappearance_rank <- NA
notincluded$torture_count <- NA
notincluded$torture_rank <- NA
notincluded$execution_count <- NA
notincluded$execution_rank <- NA
notincluded$all_commcount <- NA
notincluded$allcomm_rank <- NA
datacomp <- rbind(datacomp, notincluded)

# Mapping
## Arbitrary detention
ggplot(data = datacomp$geometry) + 
  geom_sf(aes(fill = datacomp$detention_count), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of communications \n by the UN WGAD \n (2010-2020)") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
## Disappearances
ggplot(data = datacomp$geometry) + 
  geom_sf(aes(fill = datacomp$disappearance_count), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of communications \n by the UN WGEID \n (2010-2020)") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
## Executions
ggplot(data = datacomp$geometry) + 
  geom_sf(aes(fill = datacomp$execution_count), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of communications \n by the UN SP on \n executions (2010-2020)") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
## Torture
ggplot(data = datacomp$geometry) + 
  geom_sf(aes(fill = datacomp$torture_count), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of communications \n by the UN SP on \n torture (2010-2020)") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())

## All physical integrity rights abuses
### REPRODUCES FIGURE 1 IN THE ARTICLE ###
ggplot(data = datacomp$geometry) + 
  geom_sf(aes(fill = datacomp$all_commcount), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of communications on \n physical integrity rights \n abuses (2010-2020)") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
#ggsave("Figure1.pdf", device = cairo_pdf(), width = 297, height = 210, units = "mm")

# Create Communications world map adjusted to population size 
adjmaps <- pts %>% 
  group_by(country) %>% 
  summarise(popavg = mean(population, na.rm = T), comm_count = sum(comm_count))
adjmaps$population_log <- log(adjmaps$popavg)
adjmaps$adjusted <- (adjmaps$comm_count / adjmaps$popavg)
adjmaps$adjusted_log <- (adjmaps$comm_count / adjmaps$population_log)
adjmaps$adjusted_million <- (adjmaps$adjusted * 1000000)
adjmaps <- adjmaps %>% arrange(adjusted_million)
adjmaps <- filter(adjmaps, country != "Taiwan, Province of China")
adjmaps <- filter(adjmaps, country != "Nauru")
adjmaps <- filter(adjmaps, country != "Bahrain")
adjmaps <- filter(adjmaps, country != "Maldives")
adjmaps <- filter(adjmaps, country != "Malta")
adjmaps <- filter(adjmaps, country != "Equatorial Guinea")
countries <- gisco_get_countries(year = "2016") # does only work for this year
countries <- subset(countries, CNTR_NAME != "Antarctica")
countries_small <- dplyr::select(countries, ISO3_CODE, geometry)
adjmaps$ISO3_CODE <- countrycode(adjmaps$country, origin = "country.name", destination = "iso3c")
adjmaps <- merge(adjmaps, countries_small, by = "ISO3_CODE", all.x = T, all.y = T)

# Harmonize all variables and then add all countries that are currently not included
notincluded <- dplyr::select(notincluded, ISO3_CODE, country, geometry)
notincluded$popavg <- NA
notincluded$comm_count <- NA
notincluded$population_log <- NA
notincluded$adjusted <- NA
notincluded$adjusted_log <- NA
notincluded$adjusted_million <- NA
adjmaps <- rbind(adjmaps, notincluded)

### REPRODUCES FIGURE A.1 IN THE ONLINE APPENDIX ###
ggplot(data = adjmaps$geometry) + 
  geom_sf(aes(fill = adjmaps$adjusted_million), size = 0.02) +
  scale_fill_viridis_c(option = "A", direction = -1, name = "# of communications on \n physical integrity rights \n abuses (2010-2020) \n per 1 million citizens") +
  theme(axis.text.x = element_blank(), axis.ticks = element_blank())
#ggsave("FigureA1.pdf", device = cairo_pdf(), width = 297, height = 210, units = "mm")


########################### RANK DISTANCES ###########################################

# PTS rank - Rank of all Communications added
# Positive value: PTS-rank better than ranking in # of Communications
# Negative value: PTS-rank worse than ranking in # of communications
datacomp$rankdis_all <- (datacomp$pts_rank) - (datacomp$allcomm_rank)
seeorder1 <- datacomp %>% 
  arrange(desc(rankdis_all))
seeorder1 <- seeorder1[seeorder1$country != "Gaza (Hamas)", ]
### BASED ON THIS INFORMATION TABLE A.2 & TABLE A.3 IN THE ONLINE APPENDIX HAVE BEEN CREATED ###
seeorder1 <- seeorder1[seeorder1$country != "Israel in Occupied Territories", ]
# Highest positives: United Kingdom, Australia, Switzerland, Sweden, UAE, US
# Highest negatives: Jamaica, Guinea, Cote d'Ivoire, Central African Republic, Mongolia, North Korea

### REPRODUCES FIGURE A.2 IN THE ONLINE APPENDIX ###
#pdf(file = "FigureA2.pdf")
#tiff("FigureA2.tiff", units="in", width=5, height=5, res=1200)
ggplot(aes(x = pts_rank, y = allcomm_rank), data = seeorder1) + geom_point() +  theme_bw() + geom_text(label= seeorder1$ISO3_CODE, position=position_jitter(width=4,height=4), size= 2.5) +
  xlab("Political Terror Scale rank (lower rank = more HR abuses)") +
  ylab("UNSP Communications rank (lower rank = more HR complaints)")
#dev.off()
#dev.off()


########################### PRE-ANALYSIS ###########################################

# Select relevant variables
analysisdata <- dplyr::select(pts, country, Year, comm_count, PTS, PTS_reversed, 
                              Population_log, ICCPR_rat, Armed_conflict, Elect_democracy,
                              GDPcap_log, UN_language_franca, Clean_election, Region, Liberal_democracy,
                              Alt_complaint_proc, alt_complaint_procedure_count, Judicial_ind, NHRI, CSO_repress,
                              Civil_society, OHCHR_presence, Judicial_con, y_2011, y_2012, y_2013, y_2014, y_2015,
                              y_2016, y_2017, y_2018, y_2019, y_2020)
analysisdata <- analysisdata[complete.cases(analysisdata), ]

# Descriptive statistics
### REPRODUCES TABLE A.1 IN THE ONLINE APPENDIX (only relevant covariates presented) ###
stargazer(as.data.frame(analysisdata), iqr = TRUE)

#### Check for multicollinearity ####
m1 <-
  glm.nb(comm_count ~ PTS + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election + 
           GDPcap_log + Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI,
         data = analysisdata,
         control = glm.control(maxit = 100)
  )
vif(m1)
vif_values <- vif(m1)
barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue", xlim = c(0, 10))
abline(v = 5, lwd = 3, lty = 2)


#### Check for over-dispersion ####

# Compare simple poisson with negative binomial
m1 <- glm(comm_count ~ PTS + UN_language_franca + OHCHR_presence + 
            CSO_repress + Judicial_ind + Clean_election +
            GDPcap_log + Alt_complaint_proc + Armed_conflict + 
            Population_log + NHRI,
          data = pts,
          family = "poisson"
)
summary(m1)

m2 <-
  glm.nb(comm_count ~ PTS + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           GDPcap_log + Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI,
         data = pts,
         control = glm.control(maxit = 100)
  )
summary(m2)

# Likelihood ratio test for overdispersion
L1 <- logLik(m1) # Log Likelihood of model 1
L2 <- logLik(m2) # Log Likelihood of model 2
LRT <- -2 * L1 + 2 * L2 # converges to chi^2 distribution
LRT > qchisq(0.95, df = 1) # reject null hypothesis: overdispersion present

# create variable top category of PTS
pts$pts_clean <- ifelse(pts$PTS == 1, 1, 0)



########################### ANALYSIS ###########################################


#### Negative binomial model: Standardized coefficient plot: 95% and 99% CI ####
numeric_df <- dplyr::select(pts, PTS, UN_language_franca, OHCHR_presence, 
                             GDPcap_log, CSO_repress, Judicial_ind, Clean_election,
                             Alt_complaint_proc, Armed_conflict, Population_log, NHRI)
numeric_df$country <- NULL
# Z-standardization
standardized_pts <- scale(numeric_df)
apply(standardized_pts, 2, sd, na.rm = T)
# Add dependent variable
standardized_pts <- cbind(standardized_pts, pts$comm_count)
standardized_pts <- cbind(standardized_pts, pts$year)
apply(standardized_pts, 2, sd, na.rm = T)
# As data.frame
standardized_pts <- as.data.frame(standardized_pts)
# Rename variables
standardized_pts <- rename(standardized_pts, comm_count = V12)
standardized_pts <- rename(standardized_pts, year = V13)

# Re-run models with z-standardized dataset
m1 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI,
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

m2 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m2)
s <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
stargazer(m2, se = list(s[, "Std. Error"]), type = "text")


library(broom)
library(ggplot2)
results_m1 <- m1 %>% tidy()
results_m1 <- results_m1 %>% mutate(lower = estimate - 1.96*std.error, 
                                    upper = estimate + 1.96*std.error, 
                                    lower_99 = estimate -	2.576*std.error, 
                                    upper_99 = estimate +	2.576*std.error,
                                    Scenario = "No fixed effects")

results_m2 <- m2 %>% tidy()
results_m2 <- results_m2 %>% mutate(lower = estimate - 1.96*std.error, 
                                    upper = estimate + 1.96*std.error, 
                                    lower_99 = estimate -	2.576*std.error, 
                                    upper_99 = estimate +	2.576*std.error,
                                    Scenario = "Year fixed effects")


results <- rbind(results_m1, results_m2)
results <- results %>% 
  filter(!term %in% c("factor(year)2011", "factor(year)2012", "factor(year)2013", "factor(year)2014", 
                      "factor(year)2015", "factor(year)2016", "factor(year)2017", "factor(year)2018",
                      "factor(year)2019", "factor(year)2020", "(Intercept)"))

ggplot(results) +
  geom_linerange(aes(ymin = lower_99, ymax = upper_99,
                     y = estimate, x = term, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray40", size = 0.75) +
  geom_linerange(aes(ymin = lower, ymax = upper,
                     y = estimate, x = term, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray80", size = 1.25) +
  
  geom_point(aes(y = estimate, x = term, shape = Scenario), position = position_dodge(width = 0.5), size = 2) +
  geom_hline(yintercept = 0, lty = 2) +
  xlab("") +
  coord_flip() +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  ggtitle("Predicted number of communications")


#### Negative binomial: Standardized coefficient plot with 90% and 95% CIs #####
m1 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI,
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

m2 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m2)
s <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
stargazer(m2, se = list(s[, "Std. Error"]), type = "text")


library(broom)
library(ggplot2)
results_m1 <- m1 %>% tidy()
results_m1 <- results_m1 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "No fixed effects")

results_m2 <- m2 %>% tidy()
results_m2 <- results_m2 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "Year fixed effects")


results <- rbind(results_m1, results_m2)
results <- results %>% 
  filter(!term %in% c("factor(year)2011", "factor(year)2012", "factor(year)2013", "factor(year)2014", 
                      "factor(year)2015", "factor(year)2016", "factor(year)2017", "factor(year)2018",
                      "factor(year)2019", "factor(year)2020", "(Intercept)"))

results$term[results$term == "GDPcap_log"] <- "GDP per capita"
results$term[results$term == "OHCHR_presence"] <- "OHCHR presence"
results$term[results$term == "Armed_conflict"] <- "Armed conflict"
results$term[results$term == "CSO_repress"] <- "CSO repression"
results$term[results$term == "Clean_election"] <- "Clean election"
results$term[results$term == "Alt_complaint_proc"] <- "Treaty channel"
results$term[results$term == "Judicial_ind"] <- "Jud. independence"
results$term[results$term == "UN_language_franca"] <- "UN language"
results$term[results$term == "Population_log"] <- "Population logged"
results$term[results$term == "I(PTS^2)"] <- "PTS^2"

level_order <- factor(results$term, level = c("NHRI", "Population logged", "Armed conflict", "Treaty channel", 
                                              "Clean election", "Jud. independence", "CSO repression",
                                              "GDP per capita", "UN language", "OHCHR presence", "PTS^2",
                                              "PTS"))  


### REPRODUCES FIGURE 4 IN THE ARTICLE ###
#pdf(file = "Figure4.pdf")
#tiff("Figure4.tiff", units="in", width=5, height=5, res=1200)
ggplot(results) +
  geom_linerange(aes(ymin = lower_95, ymax = upper_95,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray40", size = 0.75) +
  geom_linerange(aes(ymin = lower, ymax = upper,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray80", size = 1.25) +
  geom_point(aes(y = estimate, x = level_order, shape = Scenario), position = position_dodge(width = 0.5), size = 2) +
  geom_hline(yintercept = 0, lty = 2) +
  xlab("") +
  coord_flip() +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  ggtitle("Predicted number of UNSP Communications \n (standardized coefficients)")
#dev.off()
#dev.off()


#### Zero-inflated negative binomial: Standardized coefficient plot with 90% and 95% CIs ####
numeric_df <- dplyr::select(pts, PTS, UN_language_franca, OHCHR_presence, 
                     GDPcap_log, CSO_repress, Judicial_ind, Clean_election,
                     Alt_complaint_proc, Armed_conflict, Population_log, pts_clean, NHRI)
numeric_df$country <- NULL
# Z-standardization
standardized_pts <- scale(numeric_df)
apply(standardized_pts, 2, sd, na.rm = T)
# Add dependent variable, year, and lagged pts
standardized_pts <- cbind(standardized_pts, pts$comm_count)
standardized_pts <- cbind(standardized_pts, pts$year)
standardized_pts <- cbind(standardized_pts, pts$comm_count_lag)
apply(standardized_pts, 2, sd, na.rm = T)
# As data.frame
standardized_pts <- as.data.frame(standardized_pts)
# Rename variables
standardized_pts <- rename(standardized_pts, comm_count = V13)
standardized_pts <- rename(standardized_pts, year = V14)
standardized_pts <- rename(standardized_pts, comm_count_lag = V15)

# Run models
m1 <-
  zeroinfl(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
             GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
             Alt_complaint_proc + Armed_conflict + Population_log + NHRI | pts_clean + comm_count_lag + Armed_conflict,
           data = standardized_pts, dist = "negbin")
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

m2 <-
  zeroinfl(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
             GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
             Alt_complaint_proc + Armed_conflict + Population_log + NHRI | pts_clean + comm_count_lag + Armed_conflict +
             factor(year),
           data = standardized_pts, dist = "negbin")
summary(m2)
s <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
stargazer(m2, se = list(s[, "Std. Error"]), type = "text")

library(broom)
library(ggplot2)
library(poissonreg)
results_m1 <- m1 %>% tidy(type = "count")
results_m1 <- results_m1 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "No fixed effects")

results_m2 <- m2 %>% tidy(type = "count")
results_m2 <- results_m2 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "Year fixed effects")


results <- rbind(results_m1, results_m2)
results <- results %>% 
  filter(!term %in% c("factor(year)2011", "factor(year)2012", "factor(year)2013", "factor(year)2014", 
                      "factor(year)2015", "factor(year)2016", "factor(year)2017", "factor(year)2018",
                      "factor(year)2019", "factor(year)2020", "(Intercept)", "Log(theta)"))

results$term[results$term == "GDPcap_log"] <- "GDP per capita"
results$term[results$term == "OHCHR_presence"] <- "OHCHR presence"
results$term[results$term == "Armed_conflict"] <- "Armed conflict"
results$term[results$term == "CSO_repress"] <- "CSO repression"
results$term[results$term == "Clean_election"] <- "Clean election"
results$term[results$term == "Alt_complaint_proc"] <- "Treaty channel"
results$term[results$term == "Judicial_ind"] <- "Jud. independence"
results$term[results$term == "UN_language_franca"] <- "UN language"
results$term[results$term == "Population_log"] <- "Population logged"
results$term[results$term == "I(PTS^2)"] <- "PTS^2"

level_order <- factor(results$term, level = c("NHRI", "Population logged", "Armed conflict", "Treaty channel", 
                                              "Clean election", "Jud. independence", "CSO repression",
                                              "GDP per capita", "UN language", "OHCHR presence", "PTS^2",
                                              "PTS"))  

### REPRODUCES FIGURE A.3 IN THE ONLINE APPENDIX ###
#pdf(file = "FigureA3.pdf")
#tiff("FigureA3.tiff", units="in", width=5, height=5, res=1200)
ggplot(results) +
  geom_linerange(aes(ymin = lower_95, ymax = upper_95,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray40", size = 0.75) +
  geom_linerange(aes(ymin = lower, ymax = upper,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray80", size = 1.25) +
  geom_point(aes(y = estimate, x = level_order, shape = Scenario), position = position_dodge(width = 0.5), size = 2) +
  geom_hline(yintercept = 0, lty = 2) +
  xlab("") +
  coord_flip() +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  ggtitle("Predicted number of UNSP Communications \n (standardized coefficients)")
#dev.off()
#dev.off()



#### Negative binomial: Simulations ######

# GDP log interacted with PTS
m1 <-
  glm.nb(comm_count ~ PTS*GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election + 
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + y_2011 + y_2012 + y_2013 +
           y_2014 + y_2015 + y_2016 + y_2017 + y_2018 + y_2019 + y_2020,
         data = analysisdata,
         control = glm.control(maxit = 100)
  )
summary(m1)
screenreg(m1)

# Interaction of GDP per capita with PTS
newdata1 <- as.data.frame(cbind(intercept = 1, GDPcap_log = quantile(pts$GDPcap_log, 0.25, na.rm = T),
                                PTS = seq(min(pts$PTS, na.rm = T), max(pts$PTS, na.rm = T), length.out = 100),
                                UN_language_franca = median(pts$UN_language_franca, na.rm = T),
                                OHCHR_presence = median(pts$OHCHR_presence, na.rm = T),
                                CSO_repress = mean(pts$CSO_repress, na.rm = T),
                                Judicial_ind = mean(pts$Judicial_ind, na.rm = T),
                                Clean_election = mean(pts$Clean_election, na.rm = T),
                                Alt_complaint_proc = mean(pts$Alt_complaint_proc, na.rm = T), Armed_conflict = mean(pts$Armed_conflict),
                                Population_log = mean(pts$Population_log, na.rm = T),
                                NHRI = mean(pts$NHRI, na.rm = T),
                                y_2011 = mean(pts$y_2011), y_2012 = mean(pts$y_2012),
                                y_2013 = mean(pts$y_2013), y_2014 = mean(pts$y_2014),
                                y_2015 = mean(pts$y_2015), y_2016 = mean(pts$y_2016),
                                y_2017 = mean(pts$y_2017), y_2018 = mean(pts$y_2018),
                                y_2019 = mean(pts$y_2019), y_2020 = mean(pts$y_2020)))

predictions <- predict(m1, newdata = newdata1, type = 'response', se = T)

toplot_1<- cbind(newdata1, predictions) %>% mutate(lower = fit- 1.96*se.fit,
                                                   upper = fit + 1.96*se.fit) %>%
  mutate(scenario = "25th percentile GDP per capita (logged)")

newdata2 <- as.data.frame(cbind(intercept = 1, GDPcap_log = quantile(pts$GDPcap_log, 0.75, na.rm = T),
                                PTS = seq(min(pts$PTS, na.rm = T), max(pts$PTS, na.rm = T), length.out = 100),
                                UN_language_franca = median(pts$UN_language_franca, na.rm = T),
                                OHCHR_presence = median(pts$OHCHR_presence, na.rm = T),
                                CSO_repress = mean(pts$CSO_repress, na.rm = T),
                                Judicial_ind = mean(pts$Judicial_ind, na.rm = T),
                                Clean_election = mean(pts$Clean_election, na.rm = T),
                                Alt_complaint_proc = mean(pts$Alt_complaint_proc, na.rm = T), 
                                Armed_conflict = mean(pts$Armed_conflict),
                                Population_log = mean(pts$Population_log, na.rm = T),
                                NHRI = mean(pts$NHRI, na.rm = T),
                                y_2011 = mean(pts$y_2011), y_2012 = mean(pts$y_2012),
                                y_2013 = mean(pts$y_2013), y_2014 = mean(pts$y_2014),
                                y_2015 = mean(pts$y_2015), y_2016 = mean(pts$y_2016),
                                y_2017 = mean(pts$y_2017), y_2018 = mean(pts$y_2018),
                                y_2019 = mean(pts$y_2019), y_2020 = mean(pts$y_2020)))

predictions <- predict(m1, newdata = newdata2, type = 'response', se = T)

toplot_2 <- cbind(newdata2, predictions) %>% mutate(lower = fit- 1.96*se.fit,
                                                   upper = fit + 1.96*se.fit) %>%
  mutate(scenario = "75th percentile GDP per capita (logged)")

newdata3 <- as.data.frame(cbind(intercept = 1, GDPcap_log = quantile(pts$GDPcap_log, 0.5, na.rm = T),
                                PTS = seq(min(pts$PTS, na.rm = T), max(pts$PTS, na.rm = T), length.out = 100),
                                UN_language_franca = median(pts$UN_language_franca, na.rm = T),
                                OHCHR_presence = median(pts$OHCHR_presence, na.rm = T),
                                CSO_repress = mean(pts$CSO_repress, na.rm = T),
                                Judicial_ind = mean(pts$Judicial_ind, na.rm = T),
                                Clean_election = mean(pts$Clean_election, na.rm = T),
                                Alt_complaint_proc = mean(pts$Alt_complaint_proc, na.rm = T), 
                                Armed_conflict = mean(pts$Armed_conflict),
                                Population_log = mean(pts$Population_log, na.rm = T),
                                NHRI = mean(pts$NHRI, na.rm = T),
                                y_2011 = mean(pts$y_2011), y_2012 = mean(pts$y_2012),
                                y_2013 = mean(pts$y_2013), y_2014 = mean(pts$y_2014),
                                y_2015 = mean(pts$y_2015), y_2016 = mean(pts$y_2016),
                                y_2017 = mean(pts$y_2017), y_2018 = mean(pts$y_2018),
                                y_2019 = mean(pts$y_2019), y_2020 = mean(pts$y_2020)))

predictions <- predict(m1, newdata = newdata3, type = 'response', se = T)

toplot_3 <- cbind(newdata3, predictions) %>% mutate(lower = fit- 1.96*se.fit,
                                                    upper = fit + 1.96*se.fit) %>%
  mutate(scenario = "50th percentile GDP per capita (logged)")

toplot2 <-rbind(toplot_1, toplot_2)

### REPRODUCES FIGURE 5 IN THE ARTICLE ###
#tiff("Figure5.tiff", units="in", width=5, height=5, res=1200)
#pdf(file = "Figure5.pdf")
ggplot()+
  geom_line(data = toplot2, aes(x = PTS, y = fit, color = scenario)) +
  geom_ribbon(data = toplot2, aes(ymin=lower, ymax=upper, x=PTS, color = scenario, fill = scenario), alpha = 0.3) +
  scale_color_discrete() +
  scale_fill_discrete()+
  theme_bw() +
  labs(x= "Political Terror Scale",
       y = "Predicted Number of UNSP Communications",
       color = "",
       fill = "") +
  theme(legend.position = "bottom")
#dev.off()
#dev.off()


# OHCHR presence
m1 <-
  glm.nb(comm_count ~ PTS + GDPcap_log + OHCHR_presence + UN_language_franca +
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + y_2011 + y_2012 + y_2013 +
           y_2014 + y_2015 + y_2016 + y_2017 + y_2018 + y_2019 + y_2020,
         data = pts,
         control = glm.control(maxit = 100)
  )
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

newdata1 <- as.data.frame(cbind(intercept = 1, PTS = seq(min(pts$PTS, na.rm = T), max(pts$PTS, na.rm = T), length.out = 100),
                                GDPcap_log = mean(pts$GDPcap_log, na.rm = T),
                                OHCHR_presence = 1, UN_language_franca = median(pts$UN_language_franca, na.rm = T),
                                CSO_repress = mean(pts$CSO_repress, na.rm = T),
                                Judicial_ind = mean(pts$Judicial_ind, na.rm = T),
                                Clean_election = mean(pts$Clean_election, na.rm = T),
                                Alt_complaint_proc = mean(pts$Alt_complaint_proc, na.rm = T), Armed_conflict = mean(pts$Armed_conflict),
                                Population_log = mean(pts$Population_log, na.rm = T),
                                NHRI = mean(pts$NHRI, na.rm = T),
                                y_2011 = mean(pts$y_2011), y_2012 = mean(pts$y_2012),
                                y_2013 = mean(pts$y_2013), y_2014 = mean(pts$y_2014),
                                y_2015 = mean(pts$y_2015), y_2016 = mean(pts$y_2016),
                                y_2017 = mean(pts$y_2017), y_2018 = mean(pts$y_2018),
                                y_2019 = mean(pts$y_2019), y_2020 = mean(pts$y_2020)))

predictions <- predict(m1, newdata = newdata1, type = 'response', se = T)

toplot_1<- cbind(newdata1, predictions) %>% mutate(lower = fit- 1.96*se.fit,
                                                   upper = fit + 1.96*se.fit) %>%
  mutate(scenario = "OHCHR present")

newdata2 <- as.data.frame(cbind(intercept = 1, PTS = seq(min(pts$PTS, na.rm = T), max(pts$PTS, na.rm = T), length.out = 100),
                                GDPcap_log = mean(pts$GDPcap_log, na.rm = T),
                                OHCHR_presence = 0, UN_language_franca = median(pts$UN_language_franca, na.rm = T),
                                CSO_repress = mean(pts$CSO_repress, na.rm = T),
                                Judicial_ind = mean(pts$Judicial_ind, na.rm = T),
                                Clean_election = mean(pts$Clean_election, na.rm = T),
                                Alt_complaint_proc = mean(pts$Alt_complaint_proc, na.rm = T), Armed_conflict = mean(pts$Armed_conflict),
                                Population_log = mean(pts$Population_log, na.rm = T),
                                NHRI = mean(pts$NHRI, na.rm = T),
                                y_2011 = mean(pts$y_2011), y_2012 = mean(pts$y_2012),
                                y_2013 = mean(pts$y_2013), y_2014 = mean(pts$y_2014),
                                y_2015 = mean(pts$y_2015), y_2016 = mean(pts$y_2016),
                                y_2017 = mean(pts$y_2017), y_2018 = mean(pts$y_2018),
                                y_2019 = mean(pts$y_2019), y_2020 = mean(pts$y_2020)))

predictions2 <- predict(m1, newdata = newdata2, type = 'response', se = T)

toplot_2 <- cbind(newdata2, predictions2) %>% mutate(lower = fit- 1.96*se.fit,
                                                   upper = fit + 1.96*se.fit) %>%
  mutate(scenario = "OHCHR absent")

toplot <-rbind(toplot_2, toplot_1)

ggplot()+
  geom_line(data = toplot, aes(x = PTS, y = fit, color = scenario)) +
  geom_ribbon(data = toplot, aes(ymin=lower, ymax=upper, x=PTS, color = scenario, fill = scenario), alpha = 0.3) +
  scale_color_discrete() +
  scale_fill_discrete()+
  theme_bw() +
  labs(x= "Political Terror Scale",
       y = "Predicted Number of Communications",
       color = "",
       fill = "") +
  theme(legend.position = "bottom")


########################### REGRESSION TABLES: NEGATIVE BINOMIAL ###########################################

#### Negative binomial: Regression tables ####

m1 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s1 <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
m2 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca  + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log  + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s2 <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
m3 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           GDPcap_log + Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s3 <- coeftest(m3, vcov = vcovCL(m3, cluster = m3$country))
m4 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s4 <- coeftest(m4, vcov = vcovCL(m4, cluster = m4$country))
m5 <-
  glm.nb(comm_count ~ PTS*GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log  + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s5 <- coeftest(m5, vcov = vcovCL(m5, cluster = m5$country))

### REPRODUCES TABLE A.4 IN THE ONLNE APPENDIX ###
stargazer(list(m1, m2, m3, m4, m5), dep.var.labels = "Communications count",
          se = list(s1[, "Std. Error"], s2[, "Std. Error"], s3[, "Std. Error"],
                                              s4[, "Std. Error"], s5[, "Std. Error"]), title = "Negative binomial models",
          omit="factor\\(year", add.lines=list(c("Year fixed effects", rep("YES", 5))))


########################### ANALYSIS: LOGISTIC REGRESSION ###########################################

# Check distribution of binary communications variable
table(pts$comm_binary)

#### Logistic regression: Standardized coefficient plot with 90% and 95% CIs #####
logregvars <- dplyr::select(pts, PTS, UN_language_franca, OHCHR_presence, 
                       GDPcap_log, CSO_repress, Judicial_ind, Clean_election,
                       Alt_complaint_proc, Armed_conflict, Population_log, NHRI)
logregvars$country <- NULL
# Z-standardization
standardized_logreg <- scale(logregvars)
apply(standardized_logreg, 2, sd, na.rm = T)
# Add dependent variable
standardized_pts <- cbind(standardized_logreg, pts$comm_binary)
standardized_pts <- cbind(standardized_pts, pts$year)
apply(standardized_pts, 2, sd, na.rm = T)
# As data.frame
standardized_pts <- as.data.frame(standardized_pts)
# Rename variables
standardized_pts <- rename(standardized_pts, comm_binary = V12)
standardized_pts <- rename(standardized_pts, year = V13)

# Re-run models
m1 <- glm(comm_binary ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
            GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
            Alt_complaint_proc + Armed_conflict + 
            Population_log + NHRI, data = standardized_pts, family = binomial(link = logit))
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

m2 <- glm(comm_binary ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
            GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
            Alt_complaint_proc + Armed_conflict + 
            Population_log  + NHRI + factor(year), data = standardized_pts, family = binomial(link = logit))
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

library(broom)
library(ggplot2)
results_m1 <- m1 %>% tidy()
results_m1 <- results_m1 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "No fixed effects")

results_m2 <- m2 %>% tidy()
results_m2 <- results_m2 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "Year fixed effects")


results <- rbind(results_m1, results_m2)
results <- results %>% 
  filter(!term %in% c("factor(year)2011", "factor(year)2012", "factor(year)2013", "factor(year)2014", 
                      "factor(year)2015", "factor(year)2016", "factor(year)2017", "factor(year)2018",
                      "factor(year)2019", "factor(year)2020", "(Intercept)"))

level_order <- factor(results$term, level = c("NHRI","Population_log", "Armed_conflict", "Alt_complaint_proc", 
                                              "Clean_election", "Judicial_ind", "CSO_repress",
                                              "GDPcap_log", "UN_language_franca", "OHCHR_presence", "I(PTS^2)",
                                              "PTS"))  


ggplot(results) +
  geom_linerange(aes(ymin = lower_95, ymax = upper_95,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray40", size = 0.75) +
  geom_linerange(aes(ymin = lower, ymax = upper,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray80", size = 1.25) +
  geom_point(aes(y = estimate, x = level_order, shape = Scenario), position = position_dodge(width = 0.5), size = 2) +
  geom_hline(yintercept = 0, lty = 2) +
  xlab("") +
  coord_flip() +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  ggtitle("Predicted number of Communications \n (standardized coefficients)")


########################### REGRESSION TABLES, LOGISTIC REGRESSION ###########################################

m1 <-
  glm(comm_binary ~ PTS + I(PTS^2) + GDPcap_log + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
      family = binomial(link = logit)
  )
s1 <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
m2 <-
  glm(comm_binary ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca  + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         family = binomial(link = logit)
  )
s2 <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
m3 <-
  glm(comm_binary ~ PTS + I(PTS^2) + GDPcap_log + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           GDPcap_log + Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         family = binomial(link = logit)
  )
s3 <- coeftest(m3, vcov = vcovCL(m3, cluster = m3$country))
m4 <-
  glm(comm_binary ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         family = binomial(link = logit)
  )
s4 <- coeftest(m4, vcov = vcovCL(m4, cluster = m4$country))
m5 <-
  glm(comm_binary ~ PTS*GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         family = binomial(link = logit)
  )
s5 <- coeftest(m5, vcov = vcovCL(m5, cluster = m5$country))

### REPRODUCES TABLE A.5 IN THE ONLINE APPENDIX ###
stargazer(list(m1, m2, m3, m4, m5), se = list(s1[, "Std. Error"], s2[, "Std. Error"], s3[, "Std. Error"],
                                              s4[, "Std. Error"], s5[, "Std. Error"]), title = "Logistic regression models",
          omit="factor\\(year", add.lines=list(c("Year fixed effects", rep("YES", 5))))


#### Account for the nationality of Mandate holders ####
cor(pts$comm_count, pts$Mandate_holders) # quasi uncorrelated

# Re-run the main analysis while accounting for the nationality of mandate holders
numeric_df <- dplyr::select(pts, PTS, UN_language_franca, OHCHR_presence, 
                            GDPcap_log, CSO_repress, Judicial_ind, Clean_election,
                            Alt_complaint_proc, Armed_conflict, Population_log, NHRI, Mandate_holders)
numeric_df$country <- NULL
# Z-standardization
standardized_pts <- scale(numeric_df)
apply(standardized_pts, 2, sd, na.rm = T)
# Add dependent variable
standardized_pts <- cbind(standardized_pts, pts$comm_count)
standardized_pts <- cbind(standardized_pts, pts$year)
apply(standardized_pts, 2, sd, na.rm = T)
# As data.frame
standardized_pts <- as.data.frame(standardized_pts)
# Rename variables
standardized_pts <- rename(standardized_pts, comm_count = V13)
standardized_pts <- rename(standardized_pts, year = V14)

m1 <-
  glm.nb(comm_count ~ Mandate_holders + PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI,
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")

m2 <-
  glm.nb(comm_count ~ Mandate_holders + PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m2)
s <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
stargazer(m2, se = list(s[, "Std. Error"]), type = "text")


library(broom)
library(ggplot2)
results_m1 <- m1 %>% tidy()
results_m1 <- results_m1 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "No fixed effects")

results_m2 <- m2 %>% tidy()
results_m2 <- results_m2 %>% mutate(lower = estimate - 1.65*std.error, 
                                    upper = estimate + 1.65*std.error, 
                                    lower_95 = estimate -	1.96*std.error, 
                                    upper_95 = estimate +	1.96*std.error,
                                    Scenario = "Year fixed effects")


results <- rbind(results_m1, results_m2)
results <- results %>% 
  filter(!term %in% c("factor(year)2011", "factor(year)2012", "factor(year)2013", "factor(year)2014", 
                      "factor(year)2015", "factor(year)2016", "factor(year)2017", "factor(year)2018",
                      "factor(year)2019", "factor(year)2020", "(Intercept)"))

results$term[results$term == "GDPcap_log"] <- "GDP per capita"
results$term[results$term == "OHCHR_presence"] <- "OHCHR presence"
results$term[results$term == "Armed_conflict"] <- "Armed conflict"
results$term[results$term == "CSO_repress"] <- "CSO repression"
results$term[results$term == "Clean_election"] <- "Clean election"
results$term[results$term == "Alt_complaint_proc"] <- "Treaty channel"
results$term[results$term == "Judicial_ind"] <- "Jud. independence"
results$term[results$term == "UN_language_franca"] <- "UN language"
results$term[results$term == "Population_log"] <- "Population logged"
results$term[results$term == "I(PTS^2)"] <- "PTS^2"
results$term[results$term == "Mandate_holders"] <- "National mandate holders"


level_order <- factor(results$term, level = c("NHRI", "Population logged", "Armed conflict", "Treaty channel", 
                                              "Clean election", "Jud. independence", "CSO repression",
                                              "GDP per capita", "UN language", "OHCHR presence", "PTS^2",
                                              "PTS", "National mandate holders"))  


### REPRODUCES FIGURE A.4 IN THE APPENDIX ###
#tiff("FigureA4.tiff", units="in", width=5, height=5, res=1200)
#pdf(file = "FigureA4.pdf")
ggplot(results) +
  geom_linerange(aes(ymin = lower_95, ymax = upper_95,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray40", size = 0.75) +
  geom_linerange(aes(ymin = lower, ymax = upper,
                     y = estimate, x = level_order, shape = Scenario), 
                 position = position_dodge(width = 0.5), 
                 color = "gray80", size = 1.25) +
  geom_point(aes(y = estimate, x = level_order, shape = Scenario), position = position_dodge(width = 0.5), size = 2) +
  geom_hline(yintercept = 0, lty = 2) +
  xlab("") +
  coord_flip() +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  ggtitle("Predicted number of UNSP Communications \n (standardized coefficients)")
#dev.off()
#dev.off()

#### Robustness test with independent NHRIs ####
numeric_df <- dplyr::select(pts, PTS, UN_language_franca, OHCHR_presence, 
                            GDPcap_log, CSO_repress, Judicial_ind, Clean_election,
                            Alt_complaint_proc, Armed_conflict, Population_log, nhri_ind)
country_vec <- numeric_df$country
numeric_df$country <- NULL
# Z-standardization
standardized_pts <- scale(numeric_df)
apply(standardized_pts, 2, sd, na.rm = T)
# Add dependent variable
standardized_pts <- cbind(standardized_pts, pts$comm_count)
standardized_pts <- cbind(standardized_pts, pts$year)
apply(standardized_pts, 2, sd, na.rm = T)
# As data.frame
standardized_pts <- as.data.frame(standardized_pts)
# Rename variables
standardized_pts <- rename(standardized_pts, comm_count = V12)
standardized_pts <- rename(standardized_pts, year = V13)
m1 <-
  glm.nb(comm_count ~ PTS + I(PTS^2) + UN_language_franca + OHCHR_presence + 
           GDPcap_log + CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + nhri_ind,
         data = standardized_pts,
         control = glm.control(maxit = 100)
  )
summary(m1)
s <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
stargazer(m1, se = list(s[, "Std. Error"]), type = "text")


#### Robustness test with random effects model ####
numeric_df <- dplyr::select(pts, PTS, UN_language_franca, OHCHR_presence, 
                            GDPcap_log, CSO_repress, Judicial_ind, Clean_election,
                            Alt_complaint_proc, Armed_conflict, Population_log, NHRI)
country_vec <- numeric_df$country
numeric_df$country <- NULL
# Z-standardization
standardized_pts <- scale(numeric_df)
apply(standardized_pts, 2, sd, na.rm = T)
# Add dependent variable
standardized_pts <- cbind(standardized_pts, pts$comm_count)
standardized_pts <- cbind(standardized_pts, pts$year)
apply(standardized_pts, 2, sd, na.rm = T)
# As data.frame
standardized_pts <- as.data.frame(standardized_pts)
# Rename variables
standardized_pts <- rename(standardized_pts, comm_count = V12)
standardized_pts <- rename(standardized_pts, year = V13)
standardized_pts <- cbind(standardized_pts, country_vec)
standardized_pts$country_vec <- as.factor(standardized_pts$country_vec)
# Run random effects models
m1 <- lmer(comm_count ~ PTS + I(PTS^2) + GDPcap_log + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year) + 
         (1| country_vec), data = standardized_pts, REML = TRUE)
m2 <- lmer(comm_count ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca  + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log  + NHRI + factor(year)+ 
             (1| country_vec), data = standardized_pts, REML = TRUE)
m3 <- lmer(comm_count ~ PTS + I(PTS^2) + GDPcap_log + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           GDPcap_log + Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year)+ 
           (1| country_vec), data = standardized_pts, REML = TRUE)
m4 <- lmer(comm_count ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year) + 
           (1| country_vec), data = standardized_pts, REML = TRUE)
m5 <- lmer(comm_count ~ PTS*GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log  + NHRI + factor(year)+ 
           (1| country_vec), data = standardized_pts, REML = TRUE)
stargazer(list(m1, m2, m3, m4, m5), dep.var.labels = "Communications count",
         title = "Random effect models",
          omit="factor\\(year", add.lines=list(c("Year fixed effects", rep("YES", 5))))




########################### RE-DO ANALYSIS WITH FARISS SCORES ###########################################

# Load and prepare Fariss data
fariss_data <- read_csv('HumanRightsProtectionScores_v4.01.csv')
fariss_data <- dplyr::select(fariss_data, YEAR, COW, theta_mean)
fariss_data <- rename(fariss_data, fariss_score = theta_mean)
fariss_data <- filter(fariss_data, YEAR > 2009)
fariss_data$cow_year <- paste(fariss_data$COW, "_", fariss_data$YEAR)
fariss_data <- dplyr::select(fariss_data, cow_year, fariss_score)

# Left-join Fariss data to PTS data
pts <- left_join(pts, fariss_data, by = "cow_year")

# Re-run main models with Fariss instead of PTS
## Negative binomial regression regressions
m1 <-
  glm.nb(comm_count ~ fariss_score + I(fariss_score^2) + GDPcap_log + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log  + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s1 <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
m2 <-
  glm.nb(comm_count ~ fariss_score + I(fariss_score^2) + GDPcap_log + UN_language_franca  + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s2 <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
m3 <-
  glm.nb(comm_count ~ fariss_score + I(fariss_score^2) + GDPcap_log + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           GDPcap_log + Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s3 <- coeftest(m3, vcov = vcovCL(m3, cluster = m3$country))
m4 <-
  glm.nb(comm_count ~ fariss_score + I(fariss_score^2) + GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log  + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s4 <- coeftest(m4, vcov = vcovCL(m4, cluster = m4$country))
m5 <-
  glm.nb(comm_count ~ fariss_score*GDPcap_log + UN_language_franca + OHCHR_presence + 
           CSO_repress + Judicial_ind + Clean_election +
           Alt_complaint_proc + Armed_conflict + 
           Population_log + NHRI + factor(year),
         data = pts,
         control = glm.control(maxit = 100)
  )
s5 <- coeftest(m5, vcov = vcovCL(m5, cluster = m5$country))
stargazer(list(m1, m2, m3, m4, m5), dep.var.labels = "Communications count",
          se = list(s1[, "Std. Error"], s2[, "Std. Error"], s3[, "Std. Error"],
                    s4[, "Std. Error"], s5[, "Std. Error"]), title = "Negative binomial models",
          omit="factor\\(year", add.lines=list(c("Year fixed effects", rep("YES", 5))))


# Logistic regression models
m1 <-
  glm(comm_binary ~ fariss_score + I(fariss_score^2) + GDPcap_log + 
        CSO_repress + Judicial_ind + Clean_election +
        Alt_complaint_proc + Armed_conflict + 
        Population_log + NHRI + factor(year),
      data = pts,
      family = binomial(link = logit)
  )
s1 <- coeftest(m1, vcov = vcovCL(m1, cluster = m1$country))
m2 <-
  glm(comm_binary ~ fariss_score + I(fariss_score^2) + GDPcap_log + UN_language_franca  + 
        CSO_repress + Judicial_ind + Clean_election +
        Alt_complaint_proc + Armed_conflict + 
        Population_log + NHRI + factor(year),
      data = pts,
      family = binomial(link = logit)
  )
s2 <- coeftest(m2, vcov = vcovCL(m2, cluster = m2$country))
m3 <-
  glm(comm_binary ~ fariss_score + I(fariss_score^2) + GDPcap_log + OHCHR_presence + 
        CSO_repress + Judicial_ind + Clean_election +
        GDPcap_log + Alt_complaint_proc + Armed_conflict + 
        Population_log + NHRI + factor(year),
      data = pts,
      family = binomial(link = logit)
  )
s3 <- coeftest(m3, vcov = vcovCL(m3, cluster = m3$country))
m4 <-
  glm(comm_binary ~ fariss_score + I(fariss_score^2) + GDPcap_log + UN_language_franca + OHCHR_presence + 
        CSO_repress + Judicial_ind + Clean_election +
        Alt_complaint_proc + Armed_conflict + 
        Population_log + NHRI + factor(year),
      data = pts,
      family = binomial(link = logit)
  )
s4 <- coeftest(m4, vcov = vcovCL(m4, cluster = m4$country))
m5 <-
  glm(comm_binary ~ fariss_score*GDPcap_log + UN_language_franca + OHCHR_presence + 
        CSO_repress + Judicial_ind + Clean_election +
        Alt_complaint_proc + Armed_conflict + 
        Population_log + NHRI + factor(year),
      data = pts,
      family = binomial(link = logit)
  )
s5 <- coeftest(m5, vcov = vcovCL(m5, cluster = m5$country))
stargazer(list(m1, m2, m3, m4, m5), se = list(s1[, "Std. Error"], s2[, "Std. Error"], s3[, "Std. Error"],
                                              s4[, "Std. Error"], s5[, "Std. Error"]), title = "Logistic regression models",
          omit="factor\\(year", add.lines=list(c("Year fixed effects", rep("YES", 5))))


########################### AMELIA IMPUTATION ###########################################

library(Amelia)
set.seed(140192)
model_subdata <- dplyr::select(pts, PTS, GDPcap_log, UN_language_franca, OHCHR_presence, 
                                 CSO_repress, Judicial_ind, Clean_election,
                                 Alt_complaint_proc, Armed_conflict, comm_binary, comm_count,  
                                 Population_log, NHRI, country, Year)
model_subdata <- as.data.frame(model_subdata)
a.out <- Amelia::amelia(model_subdata, m = 5, idvars = c("comm_count", "comm_binary"), 
                        ts = "Year", cs = "country", logs = "GDPcap_log")

# Check imputations
### REPRODUCES FIGURE A.5a IN THE APPENDIX ###
#pdf(file = "FigureA5b.pdf")
#tiff("FigureA5b.tiff", units="in", width=5, height=5, res=1200)
tscsPlot(a.out, cs = "South Sudan", var = "GDPcap_log")
#dev.off()
#dev.off()
#pdf(file = "FigureA5a.pdf")
#tiff("FigureA5a.tiff", units="in", width=5, height=5, res=1200)
tscsPlot(a.out, cs = "Djibouti", var = "GDPcap_log")
#dev.off()
#dev.off()


# Listwise deletion: which cases are deleted?
model_subdata <- dplyr::select(pts, PTS, GDPcap_log, UN_language_franca, OHCHR_presence, 
                               CSO_repress, Judicial_ind, Clean_election,
                               Alt_complaint_proc, Armed_conflict, comm_binary, comm_count,  
                               Population_log, NHRI, country, Year, country_year)
used_cases <- model_subdata[complete.cases(model_subdata), ]
dropped <- dplyr::anti_join(model_subdata, used_cases, by = "country_year")
unique(dropped$country_year)
vec <- unique(dropped$country)
library(stargazer)
### REPRODUCES TABLE A.7 IN THE APPENDIX ###
stargazer(vec)
model_subdata$complete_dummy <- ifelse(complete.cases(model_subdata), 1, 0)
tapply(model_subdata$comm_count, INDEX=list(model_subdata$complete_dummy), FUN=mean)
# Higher average number of communications from complete cases

### INFORMATION USED TO CREATE TABLE A.8 IN THE APPENDIX ###
t.test(model_subdata$comm_count ~ model_subdata$complete_dummy, alternative = "two.sided")
# T-test is statistically significant


# AMELIA-imputed negative binomial models
# https://www.andrewheiss.com/blog/2018/03/07/amelia-tidy-melding/

all_imputations <- bind_rows(unclass(a.out$imputations), .id = "m") %>%
  group_by(m) %>%
  nest()
library(purrr)
library(tidyverse)

# m1
models_imputations <- all_imputations %>% mutate(model = data %>% map(~ glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + 
                                                                                 CSO_repress + Judicial_ind + Clean_election +
                                                                                 Alt_complaint_proc + Armed_conflict + 
                                                                                 Population_log + NHRI + factor(Year), data = .)),
                                                 tidied = model %>% map(~ tidy(., conf.int = TRUE)),
                                                 glance = model %>% map(~ glance(.)))

models_imputations %>%
  filter(m == "imp1") %>%
  unnest(tidied)

# Create a wide data frame of just the coefficients and standard errors
library(tidyr)
library(dplyr)
params <- models_imputations %>%
  unnest(tidied) %>% 
  dplyr::select(m, term, estimate, std.error) %>%
  gather(key, value, estimate, std.error) %>%
  spread(term, value) %>% 
  ungroup()
params 

# Extract just the coefficients
just_coefs <- params %>%
  filter(key == "estimate") %>%
  dplyr::select(-m, -key)
just_coefs

# Extract just the standard errors
just_ses <- params %>%
  filter(key == "std.error") %>%
  dplyr::select(-m, -key)
just_ses
coefs_melded <- mi.meld(just_coefs, just_ses)
coefs_melded

# Create regression summary tables
model_degree_freedom <- models_imputations %>%
  unnest(glance) %>%
  filter(m == "imp1") %>%
  pull(df.residual)

melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                      t(coefs_melded$se.mi))) %>%
  magrittr::set_colnames(c("estimate", "std.error")) %>%
  mutate(term = rownames(.)) %>%
  dplyr::select(term, everything()) %>%
  mutate(statistic = estimate / std.error,
         conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
         conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
         p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))

### REPRODUCES MODEL(1) IN TABLE A.9 IN THE APPENDIX ###
melded_summary

# m2
models_imputations <- all_imputations %>% mutate(model = data %>% map(~ glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca + 
                                                                                 CSO_repress + Judicial_ind + Clean_election +
                                                                                 Alt_complaint_proc + Armed_conflict + 
                                                                                 Population_log + NHRI + factor(Year), data = .)),
                                                 tidied = model %>% map(~ tidy(., conf.int = TRUE)),
                                                 glance = model %>% map(~ glance(.)))

models_imputations %>%
  filter(m == "imp1") %>%
  unnest(tidied)

# Create a wide data frame of just the coefficients and standard errors
library(tidyr)
params <- models_imputations %>%
  unnest(tidied) %>%
  dplyr::select(m, term, estimate, std.error) %>%
  gather(key, value, estimate, std.error) %>%
  spread(term, value) %>% 
  ungroup()
params

# Extract just the coefficients
just_coefs <- params %>%
  filter(key == "estimate") %>%
  dplyr::select(-m, -key)
just_coefs

# Extract just the standard errors
just_ses <- params %>%
  filter(key == "std.error") %>%
  dplyr::select(-m, -key)
just_ses
coefs_melded <- mi.meld(just_coefs, just_ses)
coefs_melded

# Create regression summary tables
model_degree_freedom <- models_imputations %>%
  unnest(glance) %>%
  filter(m == "imp1") %>%
  pull(df.residual)

melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                      t(coefs_melded$se.mi))) %>%
  magrittr::set_colnames(c("estimate", "std.error")) %>%
  mutate(term = rownames(.)) %>%
  dplyr::select(term, everything()) %>%
  mutate(statistic = estimate / std.error,
         conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
         conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
         p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))

### REPRODUCES MODEL(2) IN TABLE A.9 IN THE APPENDIX ###
melded_summary

# m3
models_imputations <- all_imputations %>% mutate(model = data %>% map(~ glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + OHCHR_presence + 
                                                                                 CSO_repress + Judicial_ind + Clean_election +
                                                                                 Alt_complaint_proc + Armed_conflict + 
                                                                                 Population_log + NHRI + factor(Year), data = .)),
                                                 tidied = model %>% map(~ tidy(., conf.int = TRUE)),
                                                 glance = model %>% map(~ glance(.)))

models_imputations %>%
  filter(m == "imp1") %>%
  unnest(tidied)

# Create a wide data frame of just the coefficients and standard errors
library(tidyr)
params <- models_imputations %>%
  unnest(tidied) %>%
  dplyr::select(m, term, estimate, std.error) %>%
  gather(key, value, estimate, std.error) %>%
  spread(term, value) %>% 
  ungroup()
params

# Extract just the coefficients
just_coefs <- params %>%
  filter(key == "estimate") %>%
  dplyr::select(-m, -key)
just_coefs

# Extract just the standard errors
just_ses <- params %>%
  filter(key == "std.error") %>%
  dplyr::select(-m, -key)
just_ses
coefs_melded <- mi.meld(just_coefs, just_ses)
coefs_melded

# Create regression summary tables
model_degree_freedom <- models_imputations %>%
  unnest(glance) %>%
  filter(m == "imp1") %>%
  pull(df.residual)

melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                      t(coefs_melded$se.mi))) %>%
  magrittr::set_colnames(c("estimate", "std.error")) %>%
  mutate(term = rownames(.)) %>%
  dplyr::select(term, everything()) %>%
  mutate(statistic = estimate / std.error,
         conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
         conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
         p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))

### REPRODUCES MODEL(3) IN TABLE A.9 IN THE APPENDIX ###
melded_summary


# m4
models_imputations <- all_imputations %>% mutate(model = data %>% map(~ glm.nb(comm_count ~ PTS + I(PTS^2) + GDPcap_log + UN_language_franca + OHCHR_presence + 
          CSO_repress + Judicial_ind + Clean_election +
                                        Alt_complaint_proc + Armed_conflict + 
                                        Population_log + NHRI + factor(Year), data = .)),
         tidied = model %>% map(~ tidy(., conf.int = TRUE)),
         glance = model %>% map(~ glance(.)))

models_imputations %>%
  filter(m == "imp1") %>%
  unnest(tidied)

# Create a wide data frame of just the coefficients and standard errors
library(tidyr)
params <- models_imputations %>%
  unnest(tidied) %>%
  dplyr::select(m, term, estimate, std.error) %>%
  gather(key, value, estimate, std.error) %>%
  spread(term, value) %>% 
  ungroup()
params

# Extract just the coefficients
just_coefs <- params %>%
  filter(key == "estimate") %>%
  dplyr::select(-m, -key)
just_coefs

# Extract just the standard errors
just_ses <- params %>%
  filter(key == "std.error") %>%
  dplyr::select(-m, -key)
just_ses
coefs_melded <- mi.meld(just_coefs, just_ses)
coefs_melded

# Create regression summary tables
model_degree_freedom <- models_imputations %>%
  unnest(glance) %>%
  filter(m == "imp1") %>%
  pull(df.residual)

melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                      t(coefs_melded$se.mi))) %>%
  magrittr::set_colnames(c("estimate", "std.error")) %>%
  mutate(term = rownames(.)) %>%
  dplyr::select(term, everything()) %>%
  mutate(statistic = estimate / std.error,
         conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
         conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
         p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))

### REPRODUCES MODEL(4) IN TABLE A.9 IN THE APPENDIX ###
melded_summary

# m5
models_imputations <- all_imputations %>% mutate(model = data %>% map(~ glm.nb(comm_count ~ PTS*GDPcap_log + UN_language_franca + OHCHR_presence + 
                                                                                 CSO_repress + Judicial_ind + Clean_election +
                                                                                 Alt_complaint_proc + Armed_conflict + 
                                                                                 Population_log + NHRI + factor(Year), data = .)),
                                                 tidied = model %>% map(~ tidy(., conf.int = TRUE)),
                                                 glance = model %>% map(~ glance(.)))

models_imputations %>%
  filter(m == "imp1") %>%
  unnest(tidied)

# Create a wide data frame of just the coefficients and standard errors
library(tidyr)
params <- models_imputations %>%
  unnest(tidied) %>%
  dplyr::select(m, term, estimate, std.error) %>%
  gather(key, value, estimate, std.error) %>%
  spread(term, value) %>% 
  ungroup()
params

# Extract just the coefficients
just_coefs <- params %>%
  filter(key == "estimate") %>%
  dplyr::select(-m, -key)
just_coefs

# Extract just the standard errors
just_ses <- params %>%
  filter(key == "std.error") %>%
  dplyr::select(-m, -key)
just_ses
coefs_melded <- mi.meld(just_coefs, just_ses)
coefs_melded

# Create regression summary tables
model_degree_freedom <- models_imputations %>%
  unnest(glance) %>%
  filter(m == "imp1") %>%
  pull(df.residual)

melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                      t(coefs_melded$se.mi))) %>%
  magrittr::set_colnames(c("estimate", "std.error")) %>%
  mutate(term = rownames(.)) %>%
  dplyr::select(term, everything()) %>%
  mutate(statistic = estimate / std.error,
         conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
         conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
         p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))
### REPRODUCES MODEL(5) IN TABLE A.9 IN THE APPENDIX ###
melded_summary

# END